\section{Radiance Estimation}
\label{sec:RadianceEstimation}

\subsection{Spherical Harmonics Lighting}
\label{subsec:SHLighting}
As described in \cite{Sloan:2002:PRT}, a general low frequency lighting environment can be approximated using $n$th order spherical harmonics
\begin{equation}
L(\theta,\phi)=\sum_{l=0}^{n-1}\sum_{m=-l}^{l}c_l^m Y_l^m(\theta,\phi)
\label{eqn:sh}
\end{equation}
where $Y_l^m$ are the orthogonal real spherical harmonics (SH) basis functions and $c_l^m$ are their corresponding coefficients. In this project, we use $2$nd order approximation for directional lights which requires $4$ coefficients for light representations. The volumetric radiance distribution is estimated under each lighting environment which can be represented using a SH basis function.

\subsection{Volumetric Photon Tracing}
\label{subsec:PhotonTracing}
In order to compute the global lighting for volumes, we choose volumetric photon mapping for its efficiency and simplicity. Regular mesh grids are used to store the volumetric density field where the values are defined at cell centers. Given a certain SH basis function as lighting environment (i.e. $l$ and $m$ are fixed), we just randomly emit photons from the volume boundary. For each boundary cell, a fixed number of photons are emitted with random directions $(\theta_k,\phi_k)$ with radiant energy $Y_l^m(\theta_k,\phi_k)$. This is different from the traditional volumetric photon mapping where each photon has equal amount of radiant energy. As a photon travels through the volume, we choose $\Delta\mathbf{x}=h$ as step size for photon propagation where $h$ is the grid point interval. The probability of a absorption event or scattering event happens at point $\mathbf{p}$ are integrated along the photon path from the entry point $\mathbf{p}_0$ or the location of last event
\begin{equation}
P(\mathbf{p})=1-\exp\left\{-\int_{\mathbf{p}_0}^{\mathbf{p}}(\sigma_a+\sigma_s)\rho(\mathbf{x})d\mathbf{x}\right\}
\end{equation}
where $\rho$ is density, $\sigma_a$ and $\sigma_s$ are the coefficients of absorption and scattering, respectively. If a random number $0\leq X\leq 1$ satisfies $X\leq P$, then we use Russian Roulette to determine if it is an absorption or scattering based on the probability $\sigma_a/(\sigma_a+\sigma_s)$ and $\sigma_s/(\sigma_a+\sigma_s)$. An isotropic phase function is used for scattering.
\begin{equation}
f(\theta,\phi,\theta^\prime,\phi^\prime)=\frac{1}{4\pi}
\end{equation}
The photons are cached in a balanced k-D tree at all the event locations for efficient radiance estimation.

\subsection{Radiance Estimation}
\label{subsec:RadianceEstimation}
The radiance is estimated at voxels inside the volume under each basis lighting environment. For each voxel, we collect the photons within a certain radius $r$ and evaluate the radiance
\begin{equation}
L(\theta,\phi)=\frac{\sum_k E_k W_k f(\theta,\phi,\theta_k^\prime,\phi_k^\prime)}{\frac{4}{3}\pi r^3}
\end{equation}
where $E_k$ is the photon energy, $W_k$ is the weighting kernel and $(\theta_k^\prime,\phi_k^\prime)$ is the incident angle. Currently we assume the material is isotropic therefore we only need $4$ extra variables to store the radiance at each voxel. These $4$ float-point values are then compressed and encoded using a $32$bit unsigned integer where each value is mapped from $[-1,1]\subset\mathbb{R}$ to $[0,255]\subset\mathbb{Z}$. Figure 3 shows the resultant radiance distributions from a simple low resolution test of volume photon mapping.

\begin{figure*}[ht!]
\centering
\subfigure[$l=0,m=0$]
{\includegraphics[width=0.24\linewidth]{sr00}}
\subfigure[$l=1,m=-1$]
{\includegraphics[width=0.24\linewidth]{sr10}}
\subfigure[$l=1,m=0$]
{\includegraphics[width=0.24\linewidth]{sr11}}
\subfigure[$l=1,m=1$]
{\includegraphics[width=0.24\linewidth]{sr11}}
%\vskip-.01cm
 \caption{Visualization of radiance from a certain camera view under different SH lighting basis. The radiance values are mapped to a rainbow color from blue to red.}
\label{fig:rad}
\end{figure*}
